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Deterministic Bounding Systems for Stochastic 
Compartmental Spreading Processes 

Nicholas J. Watkins, Cameron Nowzari, Victor M. Preciado, and George J. Pappas 


Abstract 

This paper studies a novel approach for approximating the behavior of compartmental spreading processes. 
In contrast to prior work, the methods developed describe a dynamics which bound the exact moment dynamics, 
without explicitly requiring a priori knowledge of non-negative (or non-positive) covariance between pairs of 
system variables. Moreover, we provide systems which provide both upper- and lower- bounds on the process 
moments. We then show that when system variables are shown to be non-negatively (or non-positively) correlated 
for all time in the system’s evolution, we may leverage the knowledge to create better approximating systems. 
We then apply the technique to several previously studied compartmental spreading processes, and compare the 
bounding systems’ performance to the standard approximations studied in prior literature. 


I. Introduction 

A. Motivation, Background, and Contributions 

Many systems of interest - viral epidemies, belief propagations, ehemieal reaetions - ean be modeled 
as eontinuous-time Markov proeesses with moments that evolve in exponentially large state spaees. 
Due to their poor scaling properties, analysis of the exact moment dynamics is prohibitive; methods of 
approximation must be invoked so as to represent the systems in a space which grows tractably with 
the size of the problem instance. Commonly, techniques referred to as moment closure approximations 
are used. 

While moment closure methods are popular and often claimed as accurate, analytic results are typically 
restricted to meager statements, e.g. demonstrating that the exact dynamics are approximated well within 
a neighborhood of the initial point. Results guaranteeing that the moments of the process remain close 
to the moments given by the approximated system often rely on simulation, and must be taken with 
caution: we begin to address this problem in our work here. 

The authors are with the Department of Electrical and Systems Engineering, University of Pennsylvania, Pennsylvania, PA 19104, USA, 
{nwatk,cnowzari,preciado,pappasg}@upenn.edu 


July 21, 2015 


DRAFT 


2 


This paper develops a method by whieh appropriate proeesses ean be approximated by a bounding system 
- one equipped with states which provably upper- and lower-bound the exact system’s moments. For 
purposes of concreteness, we focus on application to a general class of epidemic spreading processes: 
compartmental spreading models. The model we study here subsumes many of the epidemic models 
studied in literature, and so we believe to be of immediate use to the community of researchers currently 
engaged in the area. 

Literature Review: Moment closure techniques find wide use in the analysis of high-dimensional Markov 
processes, which are natural model choices in a variety of contexts. In chemistry, reactions which 
are modeled stochastically via Markov processes give rise to the Chemical Master Equations (CMEs) 
m, of which approximations are commonly studied [|2l. Similar models occur in biology, wherein 
populations are modeled stochastically [[31; in genetics, wherein protein creation and destruction are 
modeled stochastically [[H; and in epidemiology, wherein the spread of disease is modeled stochastically 

m. 

A variety of closure techniques have been used. At root, they divide into two categories: those which 
force assumptions on the process variables - such as normality (see, e.g., 0), conditional normality 
(see, e.g., [(H), and independence (see, e.g., 0, 0) - and those which construct derivative-matching 
systems near the initial system state [|3. The analytic guarantees of accuracy for such methods are not 
appropriate for problems concerning control or network design. Typical results only guarantee that the 
approximate system’s states are bounded in error over a (possibly infinitesimal) interval after the initial 
time [O. Approximation schemes exist which guarantee an error bound, such as the finite-state projection 
algorithm 0, 0, but require the repetitive solution of a system of ordinary differential equations to 
execute and only provide guarantees after the solution is computed, and so may be of limited use in 
control applications. Moreover, to the authors’ knowledge, none of these schemes provide guarantees 
as to the nature of the error: the approximate systems are not generally guaranteed to over-approximate 
or under-approximate the exact system’s moments. 

The ability to claim a rigorous upper- or lower-bound on the evolution of the process moments is 
useful for control applications. In particular, if our objective is to control an epidemic to extinction at a 
guaranteed rate, or ensure that a particular behavior survives for (or dies out within) a given amount of 
time, we must have known bounds on the moments of the process. Moreover, in order to claim control 
of an approximate model (such as done in the recent works [[91- lfmi f implies control of the modeled 
process, we must have bounding guarantees - approximations alone do not suffice. 

In compartmental spreading models, an assumption of independence is often applied (see, e.g., [(51) and 
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is referred to as a mean-field approximation. For particular choices of process structure (e.g. Susceptible- 
Infected-Susceptible SIS dill and Susceptible-Infected-Removed SIR [fT^ l. it can be shown that the 
random variables involved are non-negatively correlated IIT3ll . This supports a claim that the engendered 
approximation is an upper-bound on the infection rate of the process 0, though further claims of this 
type are lacking for other processes. In this paper, we establish a rigorous condition for testing whether 
a particular system is a bounding system for a given process, provide a construction for a non-trivial 
bounding system which is valid for any compartmental spreading process (regardless of the degree of 
correlation between the involved variables), and show that whenever knowledge of correlation between 
variables is available, it may be used to construct better bounding systems. We then apply our results to 
several models previously studied in the literature, and comment on the performance of the bounding 
systems to the approximations of prior work. 

Statement of Contributions: We define a notion of bounding systems as a technique for approximating 
stochastic processes which cannot be tractably analyzed exactly. For purposes of concreteness, we center 
our analysis on a general class of compartmental spreading models, similar to that which a general 
method of mean-field approximations was developed Q. In particular, we add the notion of “affector 
sets,” to the models studied in (Si: these allow us to cleanly consider models in which a particular 
transition is affected by neighboring nodes in a set of specified compartments. We show that for any 
model belonging to this class, a non-trivial bounding system may be constructed (Section III-AI) . and that 
whenever knowledge of the involved variables’ correlation is known a priori, it may be incorporated to 
form bounding systems which incur less error (Section fll-CI) . We demonstrate the utility of our results 
by applying them to standard epidemic models in (Section Uni). We demonstrate that whereas for some 
processes (e.g. Susceptible-Infected-Susceptible), the standard mean-field approximations are recovered 
trajectories of the bounding systems, there exist processes (e.g. Susceptible-Infected I-Susceptible- 
Infected 2- Susceptible) for which this not the case. 

B. Preliminaries 

Dynamical systems: We denote the n-dimensional Euclidean space as M". Given a vector x e M", we 
write its Euclidean norm as ||a;||. We denote a system dynamics by the state equation x = f{x), where 
/ : dom(/) !-->■ and the trajectory generated by the system as x{t), where we omit the dependence 
on time wherever it is clear from context. We say a function is globally Eipschitz continuous if there 
exists some constant L such that f{x) — f{y) ^ L\\f{x) — f{y)\\ for all x and y in the domain of /. 
The cardinality of a set Z is denoted by \Z\. 
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Probability: We denote the probability of an event A by Pr(A), the expeetation of a random variable 
X by E[X], and the indieator funetion of an event A by 1 {a}- The eovarianee of two random variables 
X and Y is defined by a(X, Y) = E[XF] - E[X]E[y]. 

Graph theory: A graph G is defined as a pair {V, E}, in whieh is a set of vertiees, and E is a set 
of edges. Two vertiees x,y e V are eonneeted if and only if {x, y) e E. 


II. COMPARTMENTAL SPREADING PROCESS MODEL 
A. Spreading Process Construction 

We eonsider a Markovian spreading proeess in whieh eaeh agent i in the spreading graph G = {V,E} 
takes membership in some eompartment c e C at all times. We will denote the membership of node i 
at time t by Xi{t) = c. Transitions from membership in eompartment c to membership in eompartment 
d (denoted by {c ^ c'}) oeeur stoehastieally via Poisson proeesses with rates where 


Yj{jeV,ceA{c,d)) l{xj(t)=c}A 


{c; c^c'} 

ij 


fc^c'} _ I Ac^c'} 
it) - i 

0 


if {c —>■ c'} is an external proeess, 
if {c —>■ c'} is an internal proeess, 
if [c —>■ d) is not possible. 


Ic- c 

and we use the notation 1{.} as the indicator of an event, ^ as the rate at which node j 

in eompartment c generates events whieh transition node i from membership in eompartment c to 
membership in eompartment d, and ^ as the rate at whieh node i generates events whieh transition 
node i from membership in eompartment c to membership in eompartment d. Note that we define a 
proeess as external whenever the transition {c —>■ c'} at node i at time t is influeneed by the eompartment 
membership of some node j at time t. A proeess is internal whenever the proeess is not external, i.e. 
when the transition at a node i is not affeeted by the eompartment membership any other node j. We 
define the set A{c, d) as the affector set assoeiated with the external transition {c —>■ d], that is the set 
of all eompartments by whieh the transition {c —>■ d] is affeeted. 

The dynamies of the first moment (i.e., the expeetation, denoted E) of this elass of proeesses ean be 
written as: 


dt 


^ (t) 

-E 

(t) 

c'£C\{c} 


- 1 

_1 


( 1 ) 


and are well defined under mild assumptions, e.g. that 0 < ^ < oo and 0 < < 00 for 

all i, j, c, c and d. We will not pause to derive these equations here, though we provide an extended 
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Fig. 1: A compartmental spreading model with \C\ = 6. Dashed arrows denote external transitions, with 
their attraetor groups eireumseribed by dashed ellipses. Solid arrows denote internal transitions. For 
this partieular example, the transitions {c —> /} and {/ d] are external, with A{c, f) = {/} and 

= {a,d}- 


diseussion of their eonstruetion in Appendix For further details, the interested reader may eonsult a 
referenee of Markov proeesses, such as lfT4l . 

This representation of the system is agent-centric: it considers states as the compartment membership 
probabilities of each node. We choose to consider this representation due to its size: it is 0{\C\n) 
dimensional. We note that this representation is not unique: we may represent the moment dynamics 
as a linear system in 0{\CY') dimensions if we consider a state space in which each possible node¬ 
compartment combination for the process is a state - this is described in much greater detail for a similar 
class of spreading models in [j5l. This approach is intractable for all but the smallest number of nodes. 
Whenever all processes are internal, the dynamics ([1]) are closed (i.e. each variable has a corresponding 
equation in the system dynamics), and may be analyzed without approximation: this representation is 
exact, which will be shown in Section ITTT-Al However, when a process {c —> c'} is external, second- 
order moments of the form E[l{ 3 .,(q=c}l{a;^(t)=c}] are introduced, preventing closure: we do not have 
equations in the system dynamics which describe the evolution of E[l{a,.(q=c}l{xj(q=c}]- Since higher 
order moments cannot, in general, be expressed in terms of first-order moments without incurring error, 
we cannot resolve this issue without approximation or a change in representation. This is primarily 
addressed in one of two ways: representing the moment dynamics as a linear system in a larger (0(|C|”)) 
state space (see [|5]| for details), or applying a moment closure technique. As the former is intractable, 
we shall concern ourselves with the latter. 
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B. Approximation Methods 

A common moment closure method in epidemic modeling is to replace the second-order moments 
E[l{xi(t)=c}l{xj(t)=c}] by a product of expectations E[l{ 2 ,,(t)=c}]IE[l{a;,■(*)=£}], in effect making an inde¬ 
pendence assumption. While for particular model choices this treatment could be a good approximation 
(see, e.g., [jH, ifTSl for discussion of the efficacy of this method as it pertains to the Susceptible-Infected- 
Susceptible model), we notice two faults which exist when applied to general models. 

Generally, this treatment can only give weak guarantees of accuracy, e.g. the existence of a neighborhood 
in which error is bounded flUl, and cannot predict whether the system over-approximates or under¬ 
approximates the model. Moreover, the systems constructed may over-approximate the model at some 
times, and under-approximate at others. This hinders a designer’s ability to relate the states of the 
approximated system to the moments of the modeled system with any confidence, preventing the 
development of rigorous guarantees with use of these approximations in many control applications, 
such as those studied in ll9l- lfTT]l . We note there are examples lHH, [fTTIl for which performing moment 
closure via an assumption of independence have resulted in approximations which estimate impossible 
probabilities, e.g. values outside of the unit interval. This, again, calls into question the efficacy of the 
technique; we will construct an alternative in what follows. 

C. Bounding Approximations 

To concentrate our analysis appropriately, we make a distinction between internal and external transitions. 
Accordingly, we let 8 (c) denote the set of all compartments c! for which either or both of the transitions 
{c ^ c'} and {c! c} are external. We let X(c) denote the set of all compartments d for which either 
or both of the transitions {c ^ d] and {d —> c} are internal. With this notation, we may rewrite the 
exact moment dynamics © as: 





dt 



E 


-E 



c'e£(c) ceA{c',c) jsV 


c'g£{c) cgA{c,c') j^V 


-t-E 


-E 



c'el{c) 


c'eX(c) 


( 2 ) 
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Note that we may evaluate the expectations and rewrite Q in terms of probabilities, which yields: 

dVi^Xi = s) 


{c; c'^c] 

ij 


dt 

c'e£{c) ceA(c',c) jeV 

+ Yj = c)d- 

c'eX(c) 


c'e£{c) ceA{c,c') jeV 


/^ /9ic; c^c 
ij 


(3) 


This representation brings us close to a system we can approximate well: we need only to bound the 
second-order moment in terms of first-order moments. For this, we recall that for any two events A and 
B, the Frechet inequalities (see,e.g., ifTSl l give that Pr(A n B) satisfies 


max{0, Pr(A) -1- Pr(i?) — 1} < Pr(v4 r^ B) ^ min{Pr(A), Pr(i?)}. 


(4) 


If we apply dH) to @ appropriately and define the state = Pr(a;j = {c}), we can verify that the 
bounds: 

dp\"^ 


dt 


^ S S 2]max{0,pf'* 

c'g£{c) cgA{c\c) jsV 
, V {c}Ac^c'} 

+ Zj Pi ^i Pi A 

c'gX(c) 


S S S ^, pf'^ }idlp 

c'e£{c) ceA(c,c') jeV 


dp, 


:c} 


dt 




c'e 6 !(c) cgA(c',c) j^V 


Z Z Z max{0, + pf ^ 

c ^ g £’( c ) csA{c^d) j^V 


,P1 .1 z.^'l _ 


c'el(c) 


(5) 


hold. We will now show how such bounds can be used to approximate the compartment membership 
probabilities of the general spreading model. 


Lemma 1 (Multivariate Comparison) Consider a system of equations 

X = f{x) 

with X E D (xMd- and f : D ^ If there exist vector functions f and f defined on D such that for 
each component fi and fi we have that 

fi{z) ^ fi{z,z), 
fi{z) > fi{z,z), 
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hold wherever Zi = Zi = Zi, Zj ^ Zj ^ Zj for all j i, and z,z,z e D. Then, the solutions to the system 

f = f{x,x), 

X = f{x,x), 

with initial conditions xifo) = x{to) = x(to), satisfy Xiit) ^ Xi{t) ^ Xiit) for all i and t ^ Iq, provided 
that the solutions x{f), x{t), and x{t) are unique and continuous. 

Proof: We argue by eontradietion. In partieular, we will show that no r e [to, oo) exists sueh that 
Xj(r) ^ Xj(r) ^ Xj(r) fails to hold. Assume for eontradietion that r' is first sueh time that the ehain of 
inequalities fails to hold. We may assume Xi(r') < xfr') for some the opposite ease will hold by a 
similar argument. By eontinuity of xft) and xff), there must exist some e > 0 sueh that Xj(r' — e) = 
Xifr' — e). Let e' be the smallest sueh value. Then, we must have that xfr' — e') < Xi{T' — e'). However, 
this eontradicts our hypothesis: at f = r' — e', it must be that x^ = ffx^x) ^ ffx) ^ fi{x,x) = Xi 
holds. Hence, no such r' exists, which implies xft) < Xi{t) ^ x, for all i and all t ^ Iq. 

Note that the case in which xfr') < xfr') for countably many i follows almost immediately. Formally, 
let X be the set of all i for which Xj(r') < Xj(r') holds. Then, by continuity of Xj(f), there exist 
Cj > 0 such that xfr' — e*) = xfr' — ef) holds. For each i, let e' be the smallest value, and define 
the index set JT” to be a permutation of X such that e'^ ^ ^ ^ Then, we must have that 

Xj^(r' — e'J < Xji(r' — e'J, but this breaks our hypothesis. The proof is completed by noting that the 
contradiction established carries through by induction. ■ 

Before continuing, it is useful to note that the functions / and / required by the hypothesis of Lemma [T] 
are not in general unique. However, it is trivially true that the pointwise minimum of all upper-bounding 
trajectories is itself an upper-bounding trajectory, and the pointwise maximum of all lower-bounding 
trajectories is itself is itself a lower-bounding trajectory: having multiple bounding systems can only 
help our analysis for any particular model. 

It is easy to see that this result gives us a mechanism for constructing a system which bounds the 
evolution of by leveraging the bound (|5l). 

Theorem 1 Given any compartmental spreading process governed by the dynamics (O, the solutions 
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of the system 


= S S S ^ “ Pi^^ ) ’ pf 

s^gS{c) cgA{c' ,c) j^V 

I /1 -{‘^}\ -{c}e{c^c^} 

+ 2j y^-Pi -Pi ^i 

c'sl(c) 


S S S max{0,pf^ + 

c'gS{c) cgA(c,c^) j^V 


p\^^ = S S S max{ 0 ,pf^ - l}/ 3 , 


c's£’(c) ceAic',c)jeV 

+ S p.'"’<*^ 

c'gX(c) 


c'g£^(c) cgA{c,c') j^V 


W}^{c'^c} _^{c}^{c^c'} 


( 6 ) 


satisfy 0 < p^^\t) ^ p'f {t) ^ ^ 1 for all t ^ to- for all p'f {Iq) = pj"^(to) = e [0,1]. 


Proof: We begin by noting that the dynamies p^ and Pi (and henee Pi) are globally Lipsehitz with 
eonstant 


L = max] 


S S Sa 

dGS{c) cgA{c'^c) JgV 


c; c - 


^c} 


+ S A 

c'eX{c) 


h'-c}| 


S S Sa 

c'gS{c) cgA{c^c') j^V 


c; c- 


>c'} 


+ S A 

c'£l(c) 




Henee, the solutions pl^\t), and pl^\t) exist, are unique, and are eontinuous. Now, we note 

pI^^ = (1 — ^ (1 “ Pi'^^) 0 and c'. It then follows that the inequalities pft) ^ 

Piit) ^ Pi{t) hold for all t ^ to by applieation of Theorem \T\ 

To show that p|'^^(t) ^ 1, we note that by eonstruetion < 0 whenever = 1, and by hypothesis 
pl'^^(to) e [0,1]. So, by eontinuity of p|‘^^(t), we have p|‘^^(t) < 1 for all t ^ to- To show that 0 < pj^^^t) 
for all t ^ to, we note that when p|^^ = 0, we have pj'^^ ^ 0. So, by p|'^^ e [0,1] and eontinuity of 

p|'^^(t), we have 0 < Pi^\t) for all t ^ to. Taken together, we have 0 < pfi) < pj(t) ^ pfi) < 1 for 

all t ^ to. ■ 


Remark 2 (Existence of Alternate Bounding Systems) Note that this eonstruetion is just one pos¬ 
sible system satisfying the hypothesis of Theorem [U In partieular, with further knowledge of the 
model’s strueture we may eonstruet alternate bounding systems, some of whieh may be more suitable 
to the problem at hand (e.g. easier to analyze or more aeeurate) than the one provided in Theorem 
[TJ As partieular examples, we demonstrate in Seetion IIII-BI that a better approximation system ean be 
eonstrueted for the Suseeptible-Infeeted-Removed (SIR) model. We show the same for the Suseeptible- 
Infeeted-Suseeptible (SIS) model in Seetion ITlI-CI Our eonstruetion testifies that one sueh system is 
possible for any eompartmental spreading proeess, but is in no means a guarantee that it is the best 
approximation. • 
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The reader may note that for all spreading process, all compartment membership probabilities at each 
node must sum to 1, and that the trajectories of dH) need not guarantee this of all possible evolutions 
allowed by the bounding trajectories. Hence, it is worthwhile to consider a means for eliminating any 
impossible state trajectories, as this may improve the approximations. 

Corollary 1 (Eliminating Impossible Trajectories) Consider the set of trajectories [p] fp} }{cec,isv} 
of dH). Then, for each c and i, the trajectory maxjO, 1 — I]c'£C\{c} under-approximator ofpf"^ 

and the trajectory min{l, 1 — '^c'eC\{c} Pi'^^} over-approximator of p’fK 

■fcl Id —Id 

Proof: This follows immediately from noting pl ^pl ^p] holds for all i and c. ■ 

We will see the effects of Corollary [T] when investigating the efficacy of the bounding system db]) in 
Section Uni 

D. Leveraging Correlations 

When considering a process for which it is known that the involved variables are non-negatively 
correlated (e.g. SIS and SIR [IT3l l. we can improve upon the upper-bounding trajectory. To show 
this formally, we will use the metric dp^^T}{f,g) = \f{t) — g{t) \ (where / and g are assumed 

continuous) as a measure for the quality of the bounding system - it tells us the size of the largest 
gap which occurs between two trajectories over the initial interval [to, ^]- Heuristically, this gives us a 
measure of the error incurred in the approximation. We can state a sufficient condition for when one 
bounding system is better than another, with respect to d{io,T}- 

Corollary 2 (Tighter Bounding Systems) Consider the metric dp^^T}{f-,g) = ^SiXte[to,T] |/(t) — gil)l 
Then, a system 

X = f{x, x) 

X = f{x, x) 
f = f{x, x) 

X = f{x, x) 

which satisfies dom(/) = dom(/) and fi{x,x) ^ fi{x,x) ^ fi{x,x) ^ fi(x,x) for all (x,x) s dom(/) 
and all {x,x) e dom(/) engenders solutions {xi{t),Xi{t),Xi{t),Xi{t)]j^i which satisfy d{to,T}{,Xi,Xi) ^ 
d{to,T}{.Xi,Xi). for all i, provided the system has unique, continuous solutions, and Xifo) = xfiR) = 
xfito) = xfto). 
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Proof: Application of Theorem [U gives us that xft) ^ xft) ^ xft) ^ Xi{t) holds for all t ^ 
and so the result follows immediately. ■ 

We may use this result to examine the effect of using knowledge of non-negative or non-positive 
eorrelation in designing a better approximating system than that whieh is eonstrueted in Theorem [B 


Theorem 3 Recall that the covariance of two random variables X and Y is defined cr(X, Y) = 
E[Xy] —E[A]E[y]. Aow, suppose that for all i, j e [n] andc,c'eC, we have a{l{xi{t)=c},i{xj{t)=c'}) ^ 
0. Then, 


Pic{xi{t) = c) FT{xj{t) = c') < min{Pr(xi(f) = c), Pr(xj(t) = c')} - (V) 

holds, and the bounding system given by 


PI _1_ _ 1 IflP’ 

ij 


Pi = Yj Y )(Pi - Y S Ximax{0,:pr- l}/3, 

dGS{c) cgA{c' ,c) j^V c'gS{c) cgA{c,c') JgV 

I /1 ~{c}\ r{c^^c} ~{c}r{c^c^} 

+ 2j Y^Pi -Pi Y 

c'el{c) 


c'g£(c) cgA(c',c) j^V c'gS{c) cgA{c^c') JgV 

{c'}r{c'^c} {c}e{c^c'} 


V P'lrP'^c} {c}.{ 

+ Zj Pi Yj - p\ 

c'£X(c) 


( 8 ) 


engenders solutions which satisfy d{to,T}iPi‘^\Pi^^) ^ d.{to,T}iPi^\Pi‘^^) far all i e \n\, c s C and all 
choices of T ^ to, where and p^f' are the solutions of the bounding system given by dH). 


Proof: The inequality © follows from the Freehet inequalities and the definition of eovarianee. 
Applieation of © shows that ® satisfies the hypothesis of Theorem |2l and so the elaim follows 
immediately. ■ 


Remark 4 (Correlations Between Some Variables) It is easy to show that a similar result holds when 
not all pairs of variables are non-negatively eorrelated: simply leave the Freehet bound substitutions 
intact wherever neeessary. A similar result holds for pairs of non-positively eorrelated variables, where 
the produet is substituted for the appropriate Freehet bound. Moreover, any a priori knowledge of non¬ 
negative (or non-positive) eorrelation allows for the construetion of an approximating system whieh is 
better with respeet to d{to,T}- • 

Theorem [3] gives us a rigorous eondition under whieh better bounding systems ean be eonstrueted. 
Whenever we have a priori knowledge of the correlation of variable pairs, we may use it to construet 
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a better system. In fact, application of the Frechet inequalities are, in effect, making an assumption 
that the correlation between variable pairs are as bad as possible, and so avoid the need for a priori 
knowledge. We may formalize this by noting that the upper bound is tight whenever the variable pair 
realizes maximal positive covariance (i.e. a{l{xi{t)=c}, '^{xj{t)=c'}) = 1 for all t), and the lower bound is 
tight whenever the variable pair realizes maximal negative covariance (i.e. '^{xj{t)=c']) = —1 

for all t). 


III. Example applications 

Having established our main technical contributions, we now demonstrate the utility of the developed 
results by applying them to several representative spreading processes. We begin by demonstrating that 
no approximation is necessary for fully internal processes (Section Illl-AI) . then continue by studying 
bounding approximations to Susceptible-Infected-Removed (SIR; Section Ull-BL Susceptible-Infected- 
Susceptible (SIS; Section ITlI-CI) . Susceptible-Infected 1-Susceptible-Infected 2-Susceptible (SRSRS; 
Section ITlI-DL and Susceptible-Exposed-Infected-Vigilant (SEIV; Section Illl-El) processes. 


A. Fully internal processes 


As a first example,we will consider a process in which all transitions are internal. We note that that the 
dynamic bounds dS]) which describe the system reduce to 


dp. 


{A 


dt 






W}.{c'^c] 


c'eX{c) 


c'gX(c) 


{c} efc^c'} 


m 


{A 


(9) 


dt 


< 




{C'}r{c'^c) 


c'el(c) 


c'gX(c) 


{c} efc^c'} 


Hence, we may consider the exact dynamics 


Id - Z 


Z Pf’«: 




( 10 ) 


c'gX(c) c'gX(c) 

directly, with no need for approximation. This demonstrates conclusively that the complexity of spreading 
processes arises from the coupling effect of external transitions - our method of approximation leads 
back to an exact representation in this case. This is a worthwhile, if obvious, result. 


B. SIR 

The Susceptible-Infected-Removed (pictured in Eigure[2l see, e.g., lfT9ll ) epidemic process has \C\ = 3, 
where the three possible compartments an agent may take are ‘Susceptible,’ (S), ‘Infected,’ (I), and ‘Re¬ 
moved,’ (R). We allow the transitions {S —> /} to be external with rate = YjjGV 
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and the transitions {/ —> R} to occur with rate We do not allow any other transitions. In 

particular, once an agent has attained state R, it will remain there ad infinitum - this allows for the 
modeling of lasting immunity after initial infection. 

It is easy to see that the bounding system dH) applies here, but we may apply Theorem [3] and the fact 
that and are non-negatively correlated [IT3l to construct the ad hoc approximating 

system: 


RS} 

Pi = 


2 max{0,pf ^ 


{I-,3^1} 


jeV 


jeV 


pp = 2(1 

jeV 


- Pi ^i 


( 11 ) 




/Vj Rj ~ -‘-id 

JeV 

pf' - pf 

which, by Corollary [2l must be as good an approximation as (jb]) with respect to and in simulation 

appears to be better. 





Fig. 2: The compartmental diagram of SIR. The dashed lines indicate external transitions; the solid 
lines indicate internal transitions. 


We note the results of a Monte Carlo simulation of the exact process as compared to the trajectories 
of the bounding system, which are given in Figure [3l From this we can see that the modeled process 
dynamics behave expected: the crude bounding system db]) gives the loosest approximation, the system 
leveraging non-negative correlation (fTTl) gives a better approximation, and the expected trajectory lies 
neatly between the over- and under-approximating trajectories of the bounding system. We note also 
that the trajectory of the original process is “tight,” to the upper and lower bounding trajectories in 
some domain for each process compartment. This is perhaps an indication that the studied bounding 
systems are close to as good as one can expect for this process, though making a formal claim to this 
end is difficult, and so this observation must be treated as a heuristic claim. In particular, it is in general 
difficult to predict a priori for which values in the domain a particular bounding trajectory will be tight. 
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(a) Expected Susceptible Fraction 



(b) Expected Infected Fraction 



(c) Expected Removed Fraction 


Fig. 3: Simulated Dynamics for the SIR Process on a 100-node random graph with connection probability 
p = 0.2. The black dashed lines are the trajectories of the bounding system ®; the red dashed lines are 
the result of applying Corollary [T] to the trajectories of ®; the magenta dashed lines are the trajectories 
of (fTT]) :the green solid line is generated from a 100-trial Monte Carlo simulation of the process. Note 
that all trajectories given are the average compartmental membership probability approximations for the 
graph. 


C. SIS 

The Susceptible-Infected-Susceptible (SIS) process (pictured in FigureS see, e.g., |[6l) epidemic process 
has |C| = 2, where the two possible compartments an agent may take are ‘Susceptible,’ (S) and ‘Infected,’ 
(I). We allow the transitions {S —> /} to be external with rate = Xijsv ^{xj{t)=i}/^ij’ 
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the transitions {/ ^ S'} to occur with rate This allows for the modeling of diseases in which 

infection may be recurrent, i.e. there is no mechanism for immunity (or removal) as in SIR. 

Some recent literature (see, e.g., IIH) has been concerned with controlling the SIS process via its mean- 
field approximation. Though an argument can be made from the combined works of dH and [|T3l that 
the mean-field approximation of SIS is an upper-bound on the first moment of the process, we can use 
our bounding system framework to study the interrelation of the bounding system method developed 
here and the standard mean-field approximation, as they relate to control of the process. 

It is easy to note that the bounding system (jb]) applies here. However, we may leverage the non-negative 
correlation of l[xi{t)=i) and for all {i,i) e V x V established in [|T3ll to write three ad hoc 

approximating systems for the SIS process: 


and 



p{I-, S^I} 
ZjjeV Pij 
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p{I} 
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Zjjsy Pij 
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= _ y rR'^ >5^ 
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= _ y rR'^ 

Z-ijeV Rij 




V aR< S^I} 

~ Z-ijeV Rij 

(i-pR)pp-p'Ar^^ 


pR) 

_ v rR' 

^jeV Rij 

max{0,y -p\''‘]-pf^A^ 

*S] 

_{S} 

Pi 

II 

1 



{S} 

Pi 

II 

1 



Rs} 

Pi 

= _ y rR’ 

Z-ijeV Rij 

max{0,pP^ - + (1 — 

AS}\Ai- 
Pi >^1 

p{S} 

= _ y rR’ 

ZjjEV Rij 

•AfA-pA) + (i-vfA 

y-5} 


1 {S} 

= i-pI 




= i-pI r 




( 12 ) 


(13) 


(14) 


Note the equations which govern the dynamics of in (fT^ and (fT3l) are identical to those of the 
standard mean-field approximation |l6l. This is an important fact. Since the mean-field trajectories are 
trajectories of a bounding system, the work done in control of SIS via mean-field approximation (see, 
e.g. a) is valid: controlling the mean-field approximation to a disease-free state must also control 
the first moment of the process to 0. However, taken together, the states of (fT^ and (fT3l) provide a 
complete picture of the system’s expected evolution, with upper and lower bounds on the probability of 
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membership of both compartments. This is an advantage of this approach over the standard mean-field 
approximation, and one which has not yet been studied in the literature. 


I 

r 

Fig. 4: The compartmental diagram of SIS. The dashed lines indicate external transitions; the solid 
lines indicate internal transitions. 

We distinguish two cases of interest in considering simulations of the process: (i) the rates are chosen 
such that an endemic (i.e. non-zero infection) steady-state exists in the mean-field approximation, and 
(ii) the rates are chosen such that the mean-field approximation dies exponentially quickly. A simulation 
of case (i) is given in Figure [H a simulation of case (ii) is given in Figure [6l 

The notable qualities of both figures are essentially the same. It is clear to see that every specified 
bounding system does, in fact, bound the evolution of the compartmental membership probabilities 
appropriately. It is interesting to note that the trajectories generated by (fT3l) appear to perform better with 
respect to d{to,T} than all other systems, although it is difficult to verify that this is the case analytically: 
the hypothesis of Corollary [2] fails. Figure [6] demonstrates that the specified bounding systems can 
perform quite well in some cases: the gap between the upper- and lower- bounding trajectories is no 
larger than 0.1 at any point in the system’s evolution. 

D. ShShS 

The Susceptible - Infected I - Susceptible - Infected 2 - Susceptible (SI 1 SI 2 S) process (pictured in 
Figure |71 see, e.g., llIOll l epidemic process has \C\ = 3, where the three possible compartments an 
agent may take are ‘Susceptible,’ (S), ‘Infected I,’ (Ji) and ‘Infected 2,’ (I 2 ). We allow the transitions 
{S Ji} to be external with rate = Yjjev the transitions {S I 2 } to 

be external with rate = Xijsy the transitions {/i ^ S'} to occur at a 

rate and the transitions {I 2 S'} to occur at a rate We allow no further transitions. 

This model allows us to consider a situation in which an agent may belong to one of two factions (/i 
or I 2 ), or be undecided S'. The standard mean-field approximation for this process can be derived by 
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Time 


(a) Expected Susceptible Fraction 
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Time 


(b) Expected Infected Fraction 


Fig. 5: Simulated Dynamics for the SIS Process on a 100-node random graph with connection probability 
p = 0.2, and rates chosen for the mean-field approximation to attain an endemic steady-state. The black 
dashed lines are the trajectories of the bounding system ® ; the magenta dashed lines are the trajectories 
of (fT^ i the blue dashed lines are the trajectories of (fT3l) : the red dashed lines are the trajectories of 
d); the green solid line is generated from a 100-trial Monte Carlo simulation of the SIS process. 
Note that all trajectories given are the average compartmental membership probability approximations 
for the graph. Note that the trajectories of (fT^ are no worse of an approximation than the trajectories 
of ® at any time: this can be verified analytically by an application of Corollary [2l 


applying results in [|5l, and is given as: 

- (1 - - 4'“') s 

jeV 

. (1 _ ( 15 ) 

jeV 

- (1 - - 4’S 

where is the mean-field approximation of the probability that node i is in compartment Ji, 
is the mean-field approximation of the probability that node i is in compartment h, and (f)] is the 
mean-field approximation of the probability that node i is in compartment S. 
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u: 0.6 
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V) 
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0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time 



(a) Expected Susceptible Fraction 



(b) Expected Infected Fraction 


Fig. 6: Simulated Dynamics for the SIS Process on a 100-node random graph with connection probability 
p = 0.2, and rates chosen for the mean-field approximation to decay exponentially. The black dashed 
lines are the trajectories of the bounding system the magenta dashed lines are the trajectories of 
(IT^ i the blue dashed lines are the trajectories of (fT3l) : the red dashed lines are the trajectories of (fT4l) : 
the green solid line is generated from a 100-trial Monte Carlo simulation of the SIS process. Note 
that all trajectories given are the average compartmental membership probability approximations for the 
graph. Note that the trajectories of (fT^ are no worse of an approximation than the trajectories of db]) 
at any time: this can be proven by application of Corollary |2] 


Some recent work IfTTl has been concerned with designing the transition rates of the network so as to 
attain a specified steady-state of the mean-field dynamics (fTSl) . however there is no work to the authors’ 
knowledge relating (fTSl) to the moment dynamics of SI 1 SI 2 S. Hence, it is interesting to consider 
simulations of the process in two cases: (i) the mean-field approximation of Ii and I 2 equilibrate to an 
endemic state, and (ii) the mean-field approximation of Ji decays exponentially to 0, while the I 2 is 
remains uncontrolled. Note that case (ii) is the control problem studied in IfTTl . 

Figure [8] studies case (i); Figure |9] studies case (ii). Note that the while the trajectories of (fTSf) do 
appear to be good approximations (heuristically), they do not systematically over- or under-estimate 
the simulated trajectory. Whether or not this is acceptable in any future practice is unclear, however 
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Fig. 7: The compartmental diagram of SI 1 SI 2 S. The dashed lines indieate external transitions; the solid 
lines indieate internal transitions. 


it does prevent a designer’s ability to make guarantees about system performanee, and as sueh it may 
be useful to eonsider the behavior of a bounding system as an alternative to mean-field approximation. 
In partieular, the eontrol of the mean-field approximation of SI 1 SI 2 S does not imply eontrol of the 
stoehastie proeess. 


E. SEIV 


The Suseeptible-Exposed-Infeeted-Vigilant (SEIV) proeess (pietured in Figure [TOl see, e.g., ifTOll l model 
has \C\ = 4, where the possible eompartments are ‘Suseeptible,’ (S), Tnfeeted,’ (I), ‘Exposed,’ (E) and 
‘Vigilant,’ (V). The transitions {S —> V}, {V S'}, {E —> /}, {/ ^ V} are speeified as internal. 
The transition {S E} is an external proeess, with transition rate ^ = Yjj '^{xj{t)=E}(3ij ’ ^ + 

A standard mean-field approximation ean be derived by applieation of results in jSll, 
and ean be written as: 

^{E}^{E^I} 

(16) 


- (1 - ^ 

jeV 


AW 


= A, 


{E^I}JE} j-{I^V},{I} 


,{V} ^ {/} ^ 


I ’ 










{v^s}AV} 


< Cl i F'\ 

where 0} ’ is an approximation of the probability that node i is a member of eompartment S, 0} ‘ is an 
approximation of the probability that node i is a member of eompartment E, is an approximation of 
the probability that node ms a member of eompartment J, and is an approximation of the probability 
that node ms a member of eompartment V. Note that (fT^ is eannot trivially be shown to be a bounding 
system; at the least, knowledge of non-negative (or non-positive) eorrelation between variables must be 
established before applying our results to derive a eomparable set of equations. However, some reeent 
work ifTOll has studied the eontrol of the SEIV proeess via the mean-field approximation (fT^ . so it is 
interesting to study the behavior of (fT^ in eomparison to a simulation of the proeess and our erude 
bounding system ®. 
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(b) Expected Infected 1 Fraction 
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(c) Expected Infected 2 Fraction 


Fig. 8: Simulated Dynamics for the SI 1 SI 2 S Process on a 100-node random graph with connection 
probability p = 0.2 with both mean-field approximations of Ji and I 2 surviving in an endemic steady 
state. The black dashed lines are the trajectories of the bounding system ® ; the red dashed lines are the 
trajectories resulting from the application of Corollary [T] to the trajectories of dH); the magenta dotted lines 
are the trajectories of the mean-field approximation (fT5l) : the green solid line is generated from a 100-trial 
Monte Carlo simulation of the process. Note that all trajectories given are the average compartmental 
membership probability approximations for the graph. Note that the mean-field approximation (fT5l) does 
not systematically over or under approximate the states. 


We again distinguish two cases of interest: (i) both the mean-field approximations and exist 
in an endemic state, and (ii) the rates are chosen such that the mean-field approximations of 0- ’ and 
exponentially decay to 0 (as was studied in ifTOl l. Case (i) is studied in Figure [HI Case (ii) is 
studied in Figure 

Here again, it is interesting to note that the mean-field dynamics ([T^ are a tighter approximation to the 
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(a) Expected Susceptible Fraction 



(b) Expected Infected 1 Fraction 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Time 


(c) Expected Infected 2 Fraction 


Fig. 9: Simulated Dynamics for the SI 1 SI 2 S Process on a 100-node random graph with connection 
probability p = 0.2 with the mean-field approximation of Ji decaying to 0 exponentially fast. The 
black dashed lines are the trajectories of the bounding system (jb]); the red dashed lines plotter are the 
trajectories resultant from applying Corollary [T] to the trajectories of ®; the magenta dotted lines are 
the trajectories of the mean-field approximation (fTSl) : the green solid line is generated from a 100-trial 
Monte Carlo simulation of the process. Note that all trajectories given are the average compartmental 
membership probability approximations for the graph. Note that the mean-field approximation (fTSl) does 
not systematically over or under approximate the states. 


simulated trajectory of the process than any of the trajectories of the bounding system, however they do 
not systematically over- or under- approximate the simulated trajectory. However, the simulations exhibit 
the property that the probabilities of membership in compartments E and / are over-approximated, 
whereas the probabilities of membership in compartments S and V are under-approximated, although 
this is a heuristic observation. This suggests that if a method for proving non-negative or non-positive 
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/ \ 



Fig. 10: The compartmental diagram of SEIV. The dashed lines indieate external transitions; the solid 
lines indicate internal transitions. 
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(a) Expected Susceptible Fraction 
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(c) Expected Infected Fraction 
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(b) Expected Exposed Fraction 



Time 


(d) Expected Vigilant Fraction 


Fig. 11: Simulated Dynamics for the SEIV Process on a 50-node random graph with connection 
probability p = 0.2. The black dashed lines are the trajectories of the bounding system ®; the red 
dashed lines are the resulting trajectories from applying Corollary [H to the trajectories of (l6l);the magenta 
dashed lines are the trajectories of the standard mean-field approximation (O; the green solid line is 
generated from a 100-trial Monte Carlo simulation of the process. Note that all trajectories given are 
the average compartmental membership probability approximations for the graph. 


correlation between process variables is developed, we may be able to construct a bounding system 
which behaves much like the mean-field approximation. 

We note that application of Corollary [T] exhibits a substantial benefit in the simulation of SEIV. 
This can be observed by examining the red dashed trajectories (those generated by the application 
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Fig. 12: Simulated Dynamics for the SEIV Process on a 50-node random graph with connection 
probability p = 0.2. The black dashed lines are the trajectories of the bounding system ®; the red 
dashed lines are the trajectories resulting from an application of Corollary [U to the trajectories of db]); 
the magenta dotted lines are the trajectories of the mean-field approximation (fT^ : the green solid line 
is generated from a 100-trial Monte Carlo simulation of the process. Note that all trajectories given are 
the average compartmental membership probability approximations for the graph. 


of Corollary [U) of the probabilities of membership in the “exposed,” compartment, which provide a 
substantial improvement over the black dashed trajectories (those generated by the dynamics (jb])). 

IV. Summary and Future Work 

Summarized most abstractly, we have shown a method constructing systems which generate trajectories 
that systematically over- and under-approximate the moments of stochastic processes with moment 
dynamics involving the expectation of a product of indicator random variables. Though widely applicable 
in the analysis of chemical, biological and engineered systems, we have focused our discussion on 
compartmental spreading processes. We have shown that a crude, but non-trivial, approximating system 
always exists for a general class of compartmental spreading process. We have shown that the crude 
construction can be improved when knowledge of non-negative or non-positive correlation between pairs 
of process variables can be shown a priori. 

In the case that all process transitions are internal, we recover the exact moment dynamics of the 
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process. In the eases of SIS and SIR, our eonstruetion leads to a reeovery of the standard mean-field 
approximation developed in [l5]| when we apply knowledge of non-negative eorrelation between infeetion 
compartment variables. In the eases of SRSRS and SEIV, we see that our erude system works as 
expeeted, but the mean-field approximations appear to agree more elosely with the simulation of the 
stoehastic proeess. 

The problem of determining whether proeess variables have non-negative or non-positive eovarianee is 
open for general eompartmental proeesses, and is of importanee. In the eontext of this work, it would 
allow us to eonstruet better bounding systems of arbitrary spreading proeesses fitting the paradigm 
presented herein without performing tedious ad hoc analyses. Though there is work near to this theme 
performed for the ease of SIS and SIR lfT3ll . it is unelear at the time of this writing that the arguments 
presented therein extend to a more general framework, sueh as the one presented in this paper. 

Though we have not diseussed it here, dynamie equations for higher order moments ean be approximated 
by the teehnique presented here as well. In the ease of eompartmental spreading processes, the seeond- 
order moment dynamies generally involve the expeetation of produets of three indieators, and we ean 
use the presented bounding teehniques almost direetly. By studying bounding systems of the first and 
second moments of a proeess, we can develop rigorous concentration arguments by applying Markov’s 
inequality, whieh would allow for a rigorous meehanism of foreeasting, however we leave further 
development to future work. 
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Appendix 

A. Construction of Exact Moment Dynamics 

From first prineiples, we derive a system of equations which describes the evolution of the probability 
that the state of node i at time t (denoted xft)) takes the value c as follows: 
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dS' [l{xi(t) = {e}}] 

dt 
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(17) 


where we have made the assumption that 

|IE[ ^ ^{xi(t+h)=c,Xi{t)=c'}] ~ ^[ ^ l{a;i(t)=c'}^j (t) h]\^hO{h), 

c'gC\{c} c'sC\{c} 

whieh holds true under mild assumptions, such as the condition 0 ^ < oo and 0 ^ /d^j’ ^ < oo 

for all i, j, c and c', which we have assumed in the body of the paper. These expressions are the 
Kolmogorov forward equations (see, e.g., [[T4l ) for the processes we are considering, where we have 
made a deliberate choice to represent the system from a microscopic (i.e. agent-centric) perspective. 
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